# deconvolution-data.jl
#
# Data matrix providing A, B, and Y.

n = 2048;  # Length of vector x
m = 256;  # Number of measurements
d = 16;  # Image dimension

# Knocked out sensor entries (these correspond to 0 in the Y vector).

knockouts = [1, 8, 17, 26, 43, 47, 49, 51, 52, 54, 56, 61, 67, 83, 84,
             95, 96, 104, 105, 106, 119, 121, 126, 136, 148, 155, 164,
             171, 173, 175, 176, 184, 192, 193, 198, 202, 204, 207,
             208, 224, 230, 238, 239, 240, 242, 244, 246, 249, 250,
             256];

# Stacked version of the observed Y matrix. All knockout entries are 0.

Y = [0.0, 1.0, 0.986736, 0.950834, 0.890618, 0.802241, 0.730228, 0.0,
     0.617109, 0.572867, 0.537012, 0.505473, 0.479585, 0.460734,
     0.434283, 0.407227, 0.0, 0.943112, 0.93028, 0.900096, 0.850174,
     0.776685, 0.715532, 0.661522, 0.616168, 0.0, 0.543681, 0.514492,
     0.490327, 0.472822, 0.446446, 0.418632, 0.917611, 0.924769,
     0.908196, 0.876034, 0.825046, 0.750903, 0.695771, 0.654471,
     0.618582, 0.585544, 0.0, 0.532771, 0.511308, 0.495924, 0.0,
     0.440252, 0.0, 0.942465, 0.0, 0.0, 0.822465, 0.0, 0.679012, 0.0,
     0.608087, 0.58723, 0.569284, 0.550538, 0.0, 0.522534, 0.496656,
     0.465714, 0.925014, 0.915113, 0.0, 0.848458, 0.79261, 0.703577,
     0.641276, 0.599525, 0.571709, 0.555844, 0.547348, 0.543127,
     0.537572, 0.534455, 0.511607, 0.479733, 0.901381, 0.889765, 0.0,
     0.0, 0.762819, 0.671031, 0.607337, 0.565367, 0.539084, 0.526613,
     0.524472, 0.529723, 0.538817, 0.54695, 0.0, 0.0, 0.870048,
     0.86508, 0.836626, 0.794755, 0.731969, 0.639567, 0.575855, 0.0,
     0.0, 0.0, 0.502627, 0.514456, 0.533332, 0.556715, 0.544042,
     0.510148, 0.88082, 0.881872, 0.852725, 0.806534, 0.734776,
     0.63639, 0.0, 0.524304, 0.0, 0.488821, 0.493562, 0.508894,
     0.533496, 0.0, 0.5615, 0.526518, 0.861447, 0.863572, 0.834812,
     0.790346, 0.719824, 0.620365, 0.551936, 0.0, 0.482313, 0.474154,
     0.481509, 0.49987, 0.528422, 0.566848, 0.573385, 0.537663,
     0.867501, 0.869613, 0.83766, 0.0, 0.71282, 0.611836, 0.54253,
     0.497902, 0.472972, 0.46596, 0.0, 0.496203, 0.527174, 0.56849,
     0.579279, 0.54319, 0.842554, 0.839698, 0.80547, 0.0, 0.686925,
     0.589376, 0.522965, 0.480936, 0.458674, 0.454589, 0.0, 0.491243,
     0.0, 0.566392, 0.0, 0.0, 0.859132, 0.850724, 0.811988, 0.766404,
     0.69285, 0.593045, 0.525301, 0.0, 0.460709, 0.45775, 0.473277,
     0.49783, 0.53043, 0.570177, 0.56848, 0.0, 0.0, 0.860428,
     0.819414, 0.77263, 0.693422, 0.0, 0.525367, 0.483245, 0.462323,
     0.0, 0.479338, 0.0, 0.534685, 0.567636, 0.0, 0.0, 0.863925,
     0.858327, 0.818639, 0.769776, 0.690254, 0.59092, 0.524347,
     0.483664, 0.464622, 0.465863, 0.487417, 0.512339, 0.537492,
     0.557952, 0.544007, 0.0, 0.884122, 0.87514, 0.830235, 0.783669,
     0.703779, 0.0, 0.534685, 0.493685, 0.47514, 0.477726, 0.50167,
     0.523405, 0.537414, 0.0, 0.0, 0.0, 0.898131, 0.0, 0.831866, 0.0,
     0.71376, 0.0, 0.542731, 0.501658, 0.0, 0.0, 0.513124, 0.526961,
     0.535028, 0.542364, 0.523565, 0.0];

# Now, let's construct the filter matrix.
#
# The translation is as follows, assuming a zero-indexed image (we
# must add one to the coordinates in Julia as it does 1-indexing):
# from coordinate (i, j) in {0, ..., d-1}^2 in the image, we get to
# vector position
#
#  k = d * i + j.
#
# Thus from vector position k, image coordinates are
#
#   j = mod(k, d).
#   i = (k - j) / d.
#
# We construct row l, corresponding to position (i, j) in the image,
# of A so that for entry k corresponding to index pair (i', j') in the
# image, we have
#
#  A[l, k] = exp(-abs(i - i') / sigma - abs(j - j') / sigma)
#
# where sigma is the bandwidth of the filter / convolution.  Note that
# Julia does 1-indexing, so we have to add 1 to l and k when we write
# A[l, k].

bandwidth = 5.0;

A = zeros(m, m);

for ll = 0:(m-1)
  j1 = mod(ll, d);
  i1 = Int((ll - j1) / d);
  for i2 = 0:(d-1)
    for j2 = 0:(d-1)
      kk = d * i2 + j2;
      A[ll+1, kk+1] = exp(-abs(i1 - i2) / bandwidth
                          - abs(j1 - j2) / bandwidth);
    end
  end
  A[ll+1, :] = A[ll+1, :] / sum(A[ll+1, :]);
end

A[A .<= 1e-5] .= 0;

A[knockouts, :] .= 0;

# Now load the B matrix.
include("B_matrix_deconv_data.jl");

